%matplotlib notebook
from ipywidgets import interact, SelectionSlider, IntSlider, FloatSlider
import matplotlib.pyplot as plt
from matplotlib import animation
import numpy as np
import scipy.signal
from IPython.display import YouTubeVideo, HTML, Audio
from bokeh.layouts import column, row
from bokeh.models import CustomJS, ColumnDataSource, Slider
from bokeh.plotting import Figure, show, output_notebook
output_notebook()
3. Efectos del muestreo y fenómeno de aliasing¶
3.1. El espectro de una señal discreta¶
3.1.1. ¿Qué le ocurre al espectro de una señal continua cuando la muestreamos?¶
Sin pérdida de generalidad, consideremos para este ejemplo la señal gaussiana
Muestrear es equivalente a multiplicar una señal continua por un tren de impulsos, también conocido como «peineta de Dirac»
donde \(F_s\) es la frecuencia de muestreo, es decir el inverso entre la separación de los dientes de la peineta
t = np.arange(-5, 5, step=0.0001)
s = lambda t, sigma=0.5 : np.exp(-0.5*t**2/sigma**2)
Fs = Slider(start=0.2, end=5, value=1, step=.01, title="Frecuencia de muestreo")
t_dirac = np.arange(-5, 5, step=1/Fs.value)
create_figure = lambda title : Figure(plot_width=280, plot_height=280, title=title,
toolbar_location="below", x_range=(-5, 5))
p1 = create_figure('Señal Gaussiana')
p2 = create_figure('Peineta de Dirac')
p3 = create_figure('Señal muestreada')
p1.line(t, s(t), line_width=3)
for p in [p1, p2, p3]:
p.xaxis[0].axis_label = 'Tiempo [s]'
source = ColumnDataSource(data=dict(t=t_dirac,
tops=np.ones_like(t_dirac),
sd=s(t_dirac)
))
p2.vbar(x='t', top='tops', source=source, width=0.05)
p3.scatter('t', 'sd', source=source)
callback = CustomJS(args=dict(source=source, Fs=Fs), code="""
var t_dirac = [];
var tops = [];
var sd = []
for (let i = -5.0; i <= 5.0; i+=1/Fs.value) {
t_dirac.push(i);
tops.push(1);
sd.push(Math.exp(-0.5*Math.pow(i, 2)/Math.pow(0.5, 2)));
}
source.data['t'] = t_dirac;
source.data['tops'] = tops;
source.data['sd'] = sd;
source.change.emit();
""")
Fs.js_on_change('value', callback)
show(column(Fs, row(p1, p2, p3)))
La transformada de Fourier del tren de impulsos es
es decir otro tren de impulsos pero en frecuencia
Pueden revisar la demostración de la transformada anterior aquí
La señal muestreada es igual a \(s(t) \cdot \upuparrows(t)\)
Para obtener el espectro de la señal muestreada recordemos que multiplicar dos señales en el dominio del tiempo corresponde a convolucionar sus transformadas de Fourier en frecuencia
Para entender la última igualdad revisemos la siguiente sección
3.1.2. ¿Qué significa convolucionar con un tren de impulsos?¶
Matemáticamente la convolución entre dos señales discretas (de una dimensión) se define como
Podemos imaginar que \(g\) se desplaza sobre \(f\) (o viceverza)
En cada paso el \(g\) desplazado se multiplica punto a punto con \(f\) y luego se suman
El resultado de convolucionar \(f\) con \(g\) es una nueva función
Si \(f\) es un tren de impulsos ocurrirá una «repetición» de \(g\) (o viceverza)
La siguiente animación muestra en la imagen superior un pulso cuadrado que se desplaza sobre un tren de impulsos en la figura superior. La figura inferior muestra el resultado de la convolución
%%capture
fig, ax = plt.subplots(2, figsize=(7, 4), sharex=True, tight_layout=True)
t = np.arange(-4, 4, step=1e-4)
def my_signal(t, a=0, T=0.5):
s = np.zeros_like(t)
s[np.absolute(t-a)<T] = 1
return s
# 2 segundos de espacio entre los dientes de la peinte
def peineta(t):
s = np.zeros_like(t)
s[::20000] = 1
return s
conv_s = np.convolve(my_signal(t), peineta(t), mode='same')
def update(a = 0):
ax[0].cla(); ax[1].cla()
p1, p2 = my_signal(t, 0.1*a - 4), peineta(t)
ax[0].plot(t, p1, label='señal');
ax[0].plot(t, p2, label='peineta');
ax[0].legend()
ax[1].plot(t, conv_s);
ax[1].set_xlabel('Tiempo [s]');
ax[1].scatter(0.1*a -4, np.sum(p1*p2), s=100, c='k')
return ()
anim = animation.FuncAnimation(fig, update, frames=80, interval=100, blit=True)
HTML(anim.to_html5_video())
Por lo tanto: Muestrear una señal en el tiempo hace que su espectro \(S(f)\), sea cual sea, se repita infinitamente. Además el espectro se repite cada \(F_s\) [Hz]
3.1.3. ¿Cúal es el espectro de la señal Gaussiana?¶
La transformada de Fourier de
es
La función gaussiana es muy especial: El espectro de una gaussiana es otra gaussiana
Puedes ver la demostración de la transformada anterior aquí
En base a la gráfica siguiente notemos que
Mientras más «ancha» sea la gaussiana en el tiempo (\(\sigma\) pequeño) más «angosto» será su espectro en frecuencia
Las gaussianas son del mismo ancho cuando \(\sigma = \frac{1}{\sqrt{2\pi}}\)
x = np.arange(-5, 5, step=0.001)
sigma = Slider(start=0.1, end=2, value=1/np.sqrt(2*np.pi), step=.01, title="sigma")
source = ColumnDataSource(data=dict(x=x,
gt=np.exp(-0.5*x**2/sigma.value**2),
gf=np.exp(-2*(np.pi*x*sigma.value)**2)*np.sqrt(2*np.pi)*sigma.value
))
create_figure = lambda title : Figure(plot_width=350, plot_height=280, title=title,
toolbar_location="below", x_range=(-5, 5))
p1 = create_figure('Dominio de tiempo', )
p2 = create_figure('Dominio de frecuencia')
p1.line('x', 'gt', source=source, line_width=3)
p2.line('x', 'gf', source=source, line_width=3)
callback = CustomJS(args=dict(source=source, sigma=sigma), code="""
var x = source.data['x'];
var gt = source.data['gt'];
var gf = source.data['gf']
for (var i = 0; i < x.length; i++) {
gt[i] = Math.exp(-0.5*Math.pow(x[i]/sigma.value, 2));
gf[i] = Math.sqrt(Math.PI*2)*Math.exp(-2*Math.pow(Math.PI*x[i]*sigma.value, 2))*sigma.value;
}
source.change.emit();
""")
sigma.js_on_change('value', callback)
show(column(sigma, row(p1, p2)))
3.1.4. ¿Cuál es el espectro de la gaussiana «discreta»?¶
Como ya vimos el espectro de una señal discreta es idéntico al de la señal continua pero «repetido» según el valor de la frecuencia de muestreo
Por ejemplo si \(\sigma=1\) y \(Fs = 2\) [Hz] tendríamos lo siguiente
def gaussiana(f, sigma):
return np.sqrt(2*np.pi*sigma**2)*np.exp(-2*(np.pi*f*sigma)**2)
def espectro_discreto(f, sigma):
S = np.zeros_like(f)
for m in range(-20, 20):
S += gaussiana(f - Fs*m, sigma)
return S
def ventana(f, Fs):
SW = np.zeros_like(f)
SW[np.absolute(f) < Fs/2] = 1
return SW
Fs = 2.
sigma = 1.
f = np.arange(-3*Fs, 3*Fs, step=1e-3)
p1 = Figure(plot_width=600, plot_height=280, toolbar_location="below")
p1.line(f, gaussiana(f, sigma), color='green', alpha=0.75,
line_width=3, legend_label='Espectro original')
p1.line(f, espectro_discreto(f, sigma), line_width=3, alpha=0.75,
legend_label='Espectro discreto')
p1.line(f, ventana(f, Fs), line_width=3, alpha=0.75, color='black',
line_dash='dashed', legend_label='Ventana cuadrada')
p1.xaxis[0].axis_label = 'Frecuencia [Hz]'
show(p1)
Por lo tanto: Si multiplicamos el espectro discreto (azul) por una ventana cuadrada de ancho \(F_s\) (negro punteado) podemos recuperar el espectro original (verde) sin pérdidas
Esta es la base del siguiente teorema fundamental
3.2. Teorema del muestreo¶
Sea una señal continua \(s(t)\) muestreada a \(F_s\) [Hz] produciendo una señal digital \(s[n] = s(t = n/F_s)\)
Si
donde
\(f_{\text{max}}\) es la componente de frecuencia más alta de la señal
\(\frac{F_s}{2}\) es la frecuencia de Nyquist
entonces la señal analógica \(s(t)\) puede ser recuperada perfectamente a partir de sus muestras discretas \(s[n]\)
Además el valor de la señal continua reconstruida es:
3.2.1. ¿Qué pasa con el espectro «discreto» si no se cumple la condición anterior?¶
Asumamos que la frecuencia de muestre se mantiene \(F_s=2\) [Hz] y que \(\sigma\) disminuye
La frecuencia máxima de la gaussiana es la «última» frecuencia donde el espectro es distinto de cero
Si \(\sigma\) disminuye la \(f_{\text{max}}\) aumenta
Fs = 2.
sigma = 0.4
f = np.arange(-3*Fs, 3*Fs, step=1e-3)
p1 = Figure(plot_width=600, plot_height=280, toolbar_location="below")
p1.line(f, gaussiana(f, sigma), color='green', alpha=0.75,
line_width=3, legend_label='Espectro original')
p1.line(f, espectro_discreto(f, sigma), line_width=3, alpha=0.75,
legend_label='Espectro discreto')
p1.line(f, ventana(f, Fs), line_width=3, alpha=0.75, color='black',
line_dash='dashed', legend_label='Ventana cuadrada')
p1.xaxis[0].axis_label = 'Frecuencia [Hz]'
show(p1)
Si asumimos que \(\sigma =1 \) se mantiene y que la frecuencia de muestreo disminuye se obtiene un efecto equivalente
Fs = 0.75
sigma = 1.
f = np.arange(-3*Fs, 3*Fs, step=1e-3)
p1 = Figure(plot_width=600, plot_height=280, toolbar_location="below")
p1.line(f, gaussiana(f, sigma), color='green', alpha=0.75,
line_width=3, legend_label='Espectro original')
p1.line(f, espectro_discreto(f, sigma), line_width=3, alpha=0.75,
legend_label='Espectro discreto')
p1.line(f, ventana(f, Fs), line_width=3, alpha=0.75, color='black',
line_dash='dashed', legend_label='Ventana cuadrada')
p1.xaxis[0].axis_label = 'Frecuencia [Hz]'
show(p1)
Debido al solapamiento en el espectro discreto (azul) se vuelve imposible recuperar el espectro original (verde) sin alteraciones. Podríamos recuperar su parte central, pero eso significa perder información de alta frecuencia
El solapamiento espectral se llama aliasing
3.3. Aliasing¶
El espectro de una señal muestreada es periódico en \(F_s\)
Si originalmente la señal tenía componentes con frecuencias mayores a \(\frac{F_s}{2}\) se produce un «traslape» o «solapamiento» espectral
Este fenomeno se llama aliasing y los componentes traslapados se denominan aliases
Cuando existe aliasing veremos que no es posible reconstruir la señal original sin ambiguedad
3.3.1. Espectro teórico de una sinusoide¶
Consideremos por ejemplo la siguiente señal sinusoidal
La transformada de Fourier de coseno es un impulso en \(f_0\) y otro en \(-f_0\)
Puedes ver la demostración de esta transformada aquí
3.3.2. Espectro de una sinusoide discreta¶
Si muestreamos a una frecuencia \(F_s\) que sea menor a \(2 f_0\) entonces habrá traslape en el espectro discreto
Por ejemplo consideremos \(f_0 = 1.23\) [Hz] y \(F_s = 2\) [Hz]. Digamos además que observamos la señal por \(100\) [s]
¿Cómo se ve el espectro de amplitud?
import scipy.fft as sfft
f0 = 1.23 # Frecuencia de la sinusoide
T = 100 # Largo temporal
Fs = 2 # Frecuencia de muestreo
ts = np.arange(0, T, step=1/Fs)
signal = lambda t, f : np.cos(2.0*np.pi*f*t)
SA = sfft.fftshift(np.absolute(sfft.fft(signal(ts, f0))))
freq = sfft.fftshift(sfft.fftfreq(n=len(ts), d=1/Fs))
p1 = Figure(plot_width=600, plot_height=280, toolbar_location="below")
p1.line(freq, SA, line_width=3)
p1.xaxis[0].axis_label = 'Frecuencia [Hz]'
show(p1)
Aparecen componentes en \(f = \pm 0.77\) [Hz]. Pero la señal era de \(1.23\) [Hz] ¿Puedes explicar porqué?
3.3.3. Espectro traslapado de la sinusoide discreta¶
Recordemos que el espectro de la señal de la señal continua se repite cada \(F_s\) [Hz] debido al muestreo
En este caso \(F_s\) es menor que \(2 f_0\) por lo que ocurrirá traslape espectral
f = np.arange(-3*Fs, 3*Fs, step=1e-3)
def espectro_coseno(f, f0):
S = np.zeros_like(f)
S[np.isclose(f, -f0)] = 1
S[np.isclose(f, f0)] = 1
return S
def espectro_discreto(f, f0):
S = np.zeros_like(f)
for m in range(-20, 20):
S += espectro_coseno(f - Fs*m, f0)
return S
def ventana(f, Fs):
SW = np.zeros_like(f)
SW[np.absolute(f) < Fs/2] = 1
return SW
p1 = Figure(plot_width=600, plot_height=280, toolbar_location="below")
p1.line(f, espectro_coseno(f, f0), color='green', alpha=0.75,
line_width=4, legend_label='Espectro original')
p1.line(f, espectro_discreto(f, f0), line_width=3, alpha=0.75,
legend_label='Espectro discreto')
p1.line(f, ventana(f, Fs), line_width=3, alpha=0.75, color='black',
line_dash='dashed', legend_label='Ventana cuadrada')
p1.xaxis[0].axis_label = 'Frecuencia [Hz]'
show(p1)
El cuadrado negro con linea punteada representa el rango
es decir lo que «podemos ver» cuando usamos la transformada de Fourier (es equivalente a la figura anterior)
3.3.4. «Aliases» de la señal¶
En este caso particular la frecuencia de la señal \(f_0\) está fuera del rango que podemos estudiar \([-F_s/2, F_s/2]\)
Además debido al solapamiento espectral aparecerán «aliases»
La frecuencia de los aliases sigue la fórmula
\(f_a = f_0 + m F_s\) [Hz] donde \(m\) es un número natural
\(f_a = m F_s - f_0\) [Hz] donde \(m\) es un número natural
Una de estas frecuencias es \(f_a = Fs - F0 = 0.773\) [Hz], que es la frecuencia que vimos en el espectro discreto anteriormente
T = 5
t = np.arange(0, T, step=1e-4)
ts = np.arange(0, T, step=1/Fs)
signal = lambda t, f : np.cos(2.0*np.pi*f*t)
p1 = Figure(plot_width=600, plot_height=280, toolbar_location="below")
p1.line(t, signal(t, f0), color='green', alpha=0.75,
line_width=3, legend_label='Señal original')
p1.scatter(ts, signal(ts, f0), size=10, legend_label='Señal muestreada')
p1.line(t, signal(t, Fs - f0), line_width=3, alpha=0.75, color='black',
line_dash='dashed', legend_label='Alias de la señal')
p1.xaxis[0].axis_label = 'Tiempo [s]'
show(p1)
El alias (linea punteada negra) se «hace pasar» por la señal original (linea verde). No podemos distinguir cual de las dos se uso para generar la señal muestreada (puntos azules). Ya no podemos reconstruir sin ambiguedad la señal
3.3.5. ¿Cómo eliminamos el aliasing?¶
Necesitamos que todas las componentes frecuenciales distintos de cero cumplan con \( |f| < \frac{Fs}{2}\)
Podemos lograr esto mediante
Filtrado: Eliminar las frecuencias mayores a \(\frac{F_s}{2}\) antes de muestrear (veremos filtrado en detalle en la unidad)
Modificar \(F_s\): Podemos aumentar \(F_s\) tal que sea dos veces mayor que la frecuencia máxima de interés, pero no siempre es fácil saber a priori cuál será dicho valor
3.4. Apéndice: Principio de incertidumbre para señales¶
Anteriormente vimos que el ancho en frecuencia y en tiempo de la gaussina son inversamentes proporcionales. Pero en realidad esto es algo que se cumple de forma general para todas las señales.
Recordemos primero el principio de incertudimbre «original», es decir el de mecámoca cuántica
El principio de incertidumbre de Heisenberg nos dice que la precisión (certeza) con que medimos la posición de una particula es inversamente proporcional a la precisión con que medimos su momentum lineal:
donde \(h\) es la constante de Planck
En señales existe un principio análogo: No podemos especificar con infinita precisión la localización temporal y frecuencial de una señal al mismo tiempo.
Denis Gabor (1946) fue el primero en darse cuenta de que el principio de incertidumbre aplica para señales.
Su teorema dice que para una señal con energía finita
con valor medio temporal
y varianza (ancho) temporal
cuya transformada de Fourier \(\mathbb{FT}[s(t)] = S(\omega)\) tiene un valor medio en frecuencia
y varianza frecuencial
se cumple que
es decir \(\Delta t\) y \(\Delta \omega\) no pueden ser arbitrariamente pequeños
El ancho temporal y el ancho frecuencial están inversamente correlacionados sin importar la señal en particular